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ABSTRACT 

Molecular  dynamics  (MD)  simulations  using  a  first 
principles-based  Embedded-atom-Method  (EAM) 
potential  are  used  to  simulate  the  exothermic  alloying 
reactions  of  a  NiAl  bilayer  at  1500  K  and  1100  K.  Both 
microcanonical  (NVE)  and  isoenthalpic  (NPH)  MD 
simulations  are  used  to  demonstrate  the  influence  of 
pressure  on  atomic  mixing  and  subsequent  alloying 
reaction.  The  NVE  simulations,  in  which  the  volume  is 
fixed  and  in  which  pressure  increases  as  the  exothermic 
reactions  cause  an  increase  in  the  system  temperature, 
have  a  much  slower  reaction  rate  than  those  of  the  NPH- 
MD  simulations  in  which  the  pressure  is  maintained  at  1 
bar.  The  mechanism  of  the  mixing  is  the  same  for  all 
simulations:  As  mixing  and  reaction  occurs  at  the 
interface,  the  heat  generated  first  melts  the  A1  layer,  and 
subsequent  mixing  leads  to  further  heat  generation  after 
which  the  Ni  layer  melts,  leading  to  additional  mixing 
until  the  alloying  reactions  are  completed.  The  results 
indicate  that  pressure  has  a  significant  influence  on  the 
rates  of  atomic  mixing  and  alloying  reactions. 

1.  INTRODUCTION 

Reactive  materials  (RMs)  hold  the  promise  of 
delivering  greater  lethality  in  future  weapons  systems  by 
serving  structural  and  energetic  functions  with  tailored 
mechanical  and  energetic  properties.  However,  to 
achieve  the  full  promise,  we  need  to  identify  and 
understand  the  fundamental  physical  and  chemical 
processes  that  control  the  energy  release  of  RMs  during 
dynamic  loading.  For  example,  we  need  to  quantify 
physical  properties  such  as  diffusion  rates  in  order  to 
guide  the  development  of  RMs  with  tailored  behavior. 
Measuring  some  of  these  properties  can  be  very 
challenging,  particularly  in  extreme  loading 
environments,  but  a  science-based  modeling  capability  to 
predict  such  properties  will  provide  a  means  for 
overcoming  measurement  limitations.  Of  particular 
importance  in  the  design  process  is  the  accurate  modeling 
of  the  dynamic  response  of  an  RM  to  rapid  thermal  or 


shock  loading,  since  they  are  the  key  components  that 
affect  the  performance  of  an  RM  in  a  weapons 
application.  One  of  the  most  widely  studied  RM  systems 
is  a  nickel-aluminum  mixture.  Studies  indicate  that  when 
such  a  system  is  used  as  a  projectile,  the  energy  release 
occurs  in  stages,  with  the  first  stage  immediately 
following  impact.  The  energy  release  is  due  to  an 
alloying  reaction  between  the  aluminum  and  nickel.  The 
second  stage  energy  release  is  the  combustion  of  the 
alloyed  metal  particles  in  the  air.  Times  scales  for  the 
various  processes  range  from  microscale  to  hundreds  of 
milliseconds.  Clearly  such  an  event  is  a  multi-scale 
phenomenon,  and  can  only  be  accurately  portrayed  using 
a  hierarchy  of  models  and  simulations  based  on 
fundamental  physical  principles.  We  are  working  toward 
instituting  a  multiscale  modeling  framework  that  can  be 
used  as  a  robust  and  reliable  modeling  and  simulation  tool 
in  the  design  and  development  process  of  RMs. 

The  most  detailed  fundamental  characterization  of 
the  chemical  and  physical  behavior  of  RMs  subjected  to 
rapid  thermal  or  shock  loading  is  obtained  through 
atomistic  simulations.  This  information  can  then  be 
transitioned  to  the  meso-  and  continuum  scale  models  that 
have  been  or  are  being  developed  for  RMs.  Such  a 
transfer  will  reduce  the  limitations  in  predictive 
capabilities  at  larger  scales  that  are  due  to  insufficient 
knowledge  of  the  details  of  the  various  chemical  and 
physical  processes  under  conditions  of  extreme 
mechanical,  thermal  and  chemical  gradients  within  a 
reacting  RM.  By  capturing  and  incorporating  the  proper 
chemical  and  physical  descriptions  of  the  dynamic 
behavior  at  the  smaller  scales  and  transitioning  this 
information  to  larger  scale  models,  more  accurate 
analyses  or  predictions  of  complex  phenomena  associated 
with  material  response  can  be  performed.  This 
information,  in  turn,  can  be  used  to  design  RMs  with 
specific  performance  for  various  weapons  that  will 
maximize  their  performance.  Incorporating  a 
fundamental  knowledge  of  material  behavior  under 
extreme  conditions  into  the  design  and  development  of 
RMs  will  substantially  reduce  the  time  and  fiscal 
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expenditures  associated  with  their  development  and 
testing. 

In  this  work,  we  present  our  efforts  to  identify  the 
role  pressure  plays  in  atomic  mixing  and  reaction  of  a 
layered  Ni-Al  system.  This  is  the  first  in  a  series  of 
studies  to  characterize  and  identify  elementary  processes 
such  as  atomic  mixing  and  melting  under  rapid  thermal 
loading,  a  phenomenon  that  will  result  from  either  rapid 
mechanical  loading  or  rapid  thermal  initiation.  We 
initiated  this  study  after  an  exploratory  microcanonical 
molecular  dynamics  (NVE-MD)  simulation  of  a  Ni/Al 
layered  system  was  performed,  in  which  the  system  was 
initially  equilibrated  to  1200K,  1  atm.  The  exothermic 
alloying  reactions  occurring  at  the  interface  produced  a 
continuous  overall  heating  and  pressure  rise  in  the 
material  to  a  point  in  time  at  which  the  alloying  reactions 
stopped.  At  the  point  of  the  quenching  of  the  reactions, 
the  system  pressure  was  extremely  high,  suggesting  that 
pressure  plays  a  significant  role  in  the  reaction  rate  for 
this  system.  Since  RMs  will  be  subjected  to  extreme 
pressures  in  military  applications,  it  is  important  to 
understand  the  dependence  on  the  atomic  mixing  and 
alloying  with  pressure. 

The  model  material  used  to  explore  these  processes  is 
a  system  similar  to  that  used  in  an  earlier  NVE-MD 
simulation  [Zhao  et  al.,  2007]  of  a  Ni/Al  bilayer  whose 
initial  temperature  was  1500  K.  In  that  study,  the 
temporal  behaviors  of  the  system  temperature  and 
pressure  were  monitored,  and  showed  very  interesting 
features,  particularly  in  the  pressure  profiles  (with  very 
similar  features  to  those  shown  in  Fig.  1(a),  in  our 
attempted  reproduction  of  the  Zhao  et  al.  results). 

The  results  of  Zhao  et  al.  [2007]  showed  that  the  system 
pressure  for  the  first  ~  75  ps  of  the  simulation  is  constant 
(~11  GPa),  after  which  it  rises  to  a  maximum  which 
coincides  with  the  time  the  Al  layer  melts.  Shortly 
thereafter,  the  Ni  layer  melts  as  heat  continues  to  evolve 
with  the  continuous  atomic  mixing/alloying.  The  pressure 
drops  and  converges  to  a  final  value  of  ~1 1.3  GPa  (in  the 
Zhao  et  al.  study)  after  alloying  and  intermixing  reaches 
completion.  In  order  to  study  the  effect  of  pressure  on  the 
overall  process,  we  will  repeat  the  same  calculation, 
except  the  equations  of  motion  will  be  coupled  with  a 
barostat,  thus  generating  a  trajectory  in  the  isoenthalpic 
ensemble  (NPH).  In  this  way,  we  can  explore  the 
temporal  behavior  of  this  system  at  a  fixed  pressure  (in 
this  case,  1  atm).  We  will  also  perform  the  simulation  at 
in  which  the  initial  bilayer  is  at  a  lower  temperature  (1100 
K),  to  examine  the  effect  of  initial  temperature  of  the 
system  on  initiation  and  propagation  of  the  alloying 
reactions. 


2.  COMPUTATIONAL  DETAILS 

The  Ni/Al  bilayer  used  in  all  simulations  is  composed 
of  47628  Ni  and  15552  Al  atoms  arranged  in  a  103  A  x  82 
A  x  89  A  filament,  with  the  Ni/Al  interface  normal  to  the 
[111]  direction,  as  described  in  Zhao  et  al.  [2007].  The 
potential  energy  function  used  in  the  simulations  is  the 
same  as  that  used  by  Zhao  et  al.  [2007]  and  was 
developed  by  Mishin  et  al.  [2002]  for  £2-NiAl.  This 
potential  follows  the  embedded-atom  method  (EAM) 
formalism  and  was  parameterized  using  both 
experimental  properties  and  ab  initio  data.  The  potential 
accurately  reproduces  the  basic  lattice  properties  of  B2- 
NiAl,  planar  faults,  and  point-defect  characteristics.  It 
also  reproduces  the  energetics  and  stability  of  all  other 
structures  included  in  the  fit.  The  FCC  lattice  constants 
used  to  generate  the  bilayer  are  3.464  A  for  Ni  and  4.0785 
A  for  Al  and  were  determined  from  NPT-MD  simulations 
of  both  bulk  Ni  and  Al  at  298  K,  1  atm. 


(b)  system  temperature  for  NVE-MD  simulation  in  which 
initial  Ni/Al  bilayer  is  at  1500  K.  The  red  line  in  Fig.  1(b) 
denotes  the  system  temperature  from  NPH-MD 

simulations  under  the  same  initial  conditions.  The  inset  in 
Fig.  1(b)  is  an  enlargement  of  these  system  temperature 
profiles  between  0  and  60  ps.  The  green  line  in  Fig.  1(b) 
denotes  the  system  temperature  from  NPH-MD 

simulations  for  a  bilayer  with  initial  temperature  of  1 1 00 
K. 
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Few  details  of  the  initial  conditions  used  by  Zhao  et 
al  in  their  NVE-MD  simulations  were  given,  other  than 
the  system  was  assigned  an  initial  temperature  of  1500  K. 
Thus,  we  arbitrarily  chose  to  initiate  the  trajectory  from 
an  unrelaxed  model  in  which  all  atoms  are  in  the  FCC 
lattice  sites.  Atomic  velocities  were  randomly  assigned  to 
correspond  to  an  initial  temperature  of  1500  K,  and  a 
short  NVE-MD  trajectory  was  run  in  which  the  atomic 
velocities  were  scaled  at  every  0.0072  ps  to  maintain  an 
instantaneous  system  temperature  of  1500  K.  We  were 
unable  to  perform  an  extended  equilibration  without 
alloying  reaction  at  the  interface  occurring;  thus,  after  a 
short  period  (after  which  the  system  was  not  completely 
equilibrated),  we  began  a  “production”  NVE-MD 
simulation  for  a  duration  450  ps  to  compare  with  the  Zhao 
et  al.  results  [see  Fig.  1(a)]. 

Initial  conditions  for  the  NPH  simulations  were 
generated  after  a  series  of  NPT-MD  equilibration 
trajectories  to  raise  the  system  temperature  gradually. 
First,  the  same  Ni/Al  bilayer  used  for  the  previously 
described  NVE-MD  simulation  was  equilibrated  using 
NPT-MD  simulations  at  T=298  K,  P  =  1  atm.  Upon 
completion  of  this  trajectory,  the  final  state  was  used  in  a 
subsequent  short  NPT-MD  simulation  (15  ps,  1  atm)  with 
a  target  temperature  of  T=800  K.  No  reaction  occurred  at 
the  interface  by  the  end  of  this  second  trajectory.  This 
procedure  was  repeated  choosing  a  target  temperature  of 
1100  K,  and  the  system  configuration  at  the  end  of  the 
equilibration  trajectory  was  used  for  initial  conditions  of 
the  “production”  NPH-MD  simulation  at  T=1100K. 
Again,  reactions  at  the  interface  had  not  yet  begun  at  the 
end  of  the  equilibration  trajectory  at  1100  K.  We 
attempted  to  ramp  the  temperature  up  to  1500  K  using  the 
aforementioned  procedure,  and  equilibrate  to  P=latm; 
however,  reaction  at  the  interface  began  almost 
immediately.  Thus,  in  order  to  study  the  system  at  1500 
K,  1  atm,  we  were  forced  to  use  an  incompletely 
equilibrated  system  configuration  for  which  only  a  few 
equilibration  steps  were  taken  to  provide  the  proper 
system  temperature  and  an  approximately  relaxed  system 
pressure.  Although  not  fully  equilibrated  this  initial  state 
was  sufficiently  relaxed  that  the  imposed  pressure  (1  atm) 
was  quickly  reached. 


3.  RESULTS  AND  DISCUSSION 

3.1  NVE-MD,  T=1500  K.  The  behavior  of  the  pressure 
as  a  function  of  time  during  this  trajectory  is  shown  in 
Fig.  1(a)  and  can  be  directly  compared  with  Fig.  5  of 
Zhao  et  al.  [2007].  The  corresponding  system 
temperature  is  shown  in  Fig.  1(b).  The  features  of  the 
pressure  profile  are  very  similar  to  those  of  Zhao  et  al. 
[2007],  except  they  are  shifted  in  time  by  ~  20  ps,  and  the 
pressure  profile  is  consistently  higher  by  0.3  GPa  for  the 
duration  of  the  trajectory.  Zhao  et  al.  [2007]  did  not 


provide  a  system  temperature  profile  in  their  study; 
however,  they  provided  a  figure  that  described  the 
temporal  behavior  of  the  Al  layer  only.  Unfortunately, 
Zhao  et  al.  [2007]  did  not  define  what  the  Al  layer  is,  and 
we  were  unable  to  determine  whether  it  is  a  fixed  spatial 
range  or  whether  it  was  the  portion  of  the  material 
composed  only  of  Al  atoms.  Regardless,  we  were  unable 
to  reproduce  that  figure.  In  Fig.  1(b),  however,  we  see 
that  the  system  temperature  begins  to  increase  at  about  75 
ps,  and  there  is  a  very  slight  dip  at  -160  ps,  which 
corresponds  to  the  peak  pressure  of  the  system  during  this 
trajectory  (Fig.  la).  This  is  the  point  at  which  Al  melts, 
and  corresponds  to  a  system  temperature  of  -  1800  K. 
After  a  short  time  and  small  temperature  increase,  the 
temperature  profile  shows  yet  another  “dip”,  which 
corresponds  to  the  point  at  which  Ni  melts  in  this 
simulation. 

3.2  NPH-MD,  T=1500  K.  The  temporal  behavior  of  the 
system  temperature  in  this  trajectory  is  also  shown  in  Fig. 
1(b)  for  comparison  with  the  results  of  the  NVE-MD 
trajectory.  In  the  NPH-MD  simulation,  pressure  is 
continually  relieved  as  the  system  temperature  rises  due  to 
the  coupling  of  the  equations  of  motion  with  a  barostat 
that  requires  the  system  to  maintain  an  imposed  pressure 
(in  this  case,  1  atm).  The  profiles  are  dramatically 
different.  The  temperature  of  the  NPH-MD  simulation 
immediately  drops  during  the  first  few  ps  of  the  trajectory 
integration  (see  inset),  which  corresponds  to  the  almost 
immediate  melting  of  the  Al  layer  at  a  temperature  of  - 
1500  K.  The  system  temperature  then  steadily  increases 
until  -  30  ps,  at  which  time  there  is  another  slight  dip  in 
temperature  before  the  temperature  continually  rises  and 
levels  off  to  2164  K,  -  100  K  lower  than  that  of  the 
NVE-MD  simulations.  This  second  dip  corresponds  to 
the  melting  of  the  Ni  layer,  at  a  temperature  of  -  1600  K. 
Both  melting  temperatures  are  substantially  lower  than 
those  in  the  system  constrained  to  a  constant  volume  (i.e. 
high  pressure).  The  NPH-MD  temperature  profiles 
indicate  that  the  alloying  reactions  and  mixing  have 
reached  completion  by  200  ps,  well  before  that  of  the 
NVE-MD  simulations. 

A  comparison  of  the  snapshots  from  the  NPH-MD 
and  NVE-MD  trajectories  show  no  visual  differences  in 
the  materials  during  the  transition  from  bilayer  to  fully- 
mixed  material;  it  appears  the  most  significant  impact  of 
the  pressure  is  on  the  rates  of  alloying,  mixing  and 
melting. 

3.3  NPH-MD,  T=1100  K.  NPH-MD  simulations  at  this 
lower  temperature  produced  a  different  system 
temperature  profile  than  the  NPH-MD  simulations  at 
T=1500  K,  as  seen  in  Fig.  1(b),  and  provided  a  profile  that 
was  similar  to  the  NVE-MD  simulations  at  1500  K.  In 
this  trajectory,  the  system  temperature  fluctuates  about 
1100  K  up  to  100  ps,  after  which  time  the  temperature 
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steadily  rises  to  -1200  K  at  190  ps  [see  inset  of  Fig.  2(a)], 
which  is  the  point  at  which  the  A1  layer  melts.  A  small 
dip  in  the  system  temperature  occurs  at  -  275  ps,  the  point 
at  which  the  Ni  layer  melts  (~  1480  K).  The  temperature 
then  continually  increases  for  the  next  150  ps  before 
leveling  off  to  a  final  value  of  1764  K. 

Although  we  are  showing  system  temperature  profiles  for 
the  three  trajectories,  the  material  does  not  necessarily 
have  a  uniform  temperature  distribution  throughout  the 
material.  To  explore  the  temperature  distribution  in  the 
material  as  it  undergoes  alloying,  we  will  examine  first 
the  temporal  profiles  of  layers  of  the  material  during  the 
NPH-MD  simulation  in  which  the  initial  temperature  is 
1100  K.  We  partitioned  the  bilayer  into  12  sections,  for 
which  temperature  within  each  section  was  determined. 
Figure  2  shows  the  temporal  profiles  for  each  section  of 
the  material  for  the  duration  of  the  trajectory. 


Figure  2.  Temperature  as  a  function  of  time  at  partitions 
of  the  layered  material.  The  left-most  schematic  depicts 
the  overall  Ni/Al  bilayer  and  the  12  partitions  within 
which  temperature  of  the  material  was  averaged.  The 
corresponding  temperature  profiles  of  the  material  are  at 
the  right  of  the  schematic. 

For  convenience,  the  left-most  shaded  area  in  Fig.  2 
is  a  schematic  of  the  initial  Al/Ni  bilayer,  illustrating  the 
partitioning  of  the  system  to  explore  local  temperature 
throughout  the  trajectory.  For  the  most  part,  the  overall 
features  of  the  layers  track  those  of  the  system 
temperature,  with  small  differences  corresponding  to  the 


features  associated  with  melting.  For  example,  layers  3 
and  9  correspond  to  the  central  regions  of  the  A1  and  Ni 
layers,  respectively;  these  show  the  sharpest  “dips”  in  the 
temperature  profiles  that  are  associated  with  their 
respective  melting  points. 

Snapshots  of  the  bilayer  at  232,  275  and  710  ps 
during  this  trajectory  are  shown  in  Fig.  3.  Under  each 
snapshot  is  the  corresponding  temperature  contour  of  the 
material  relative  to  the  instantaneous  system  temperature. 
It  is  clear  that  at  232  ps,  the  aluminum  layer  is  disordered, 
and  the  melted  A1  layer  is  much  hotter  than  the  Ni  layer, 
indicating  that  mixing  into  the  melted  A1  is  occurring 
rapidly. 


Figure  3.  Snapshots  at  232  ps 
(right),  275  ps  (middle)  and  710  ps 
(left).  Corresponding  temperature 
contours  of  the  system  relative  to 
the  system  temperature  are 
underneath  each  snapshot. 


At  275  ps,  the  Ni  region  of  the  bilayer  has  a  semblance  of 
order  but  the  temperature  contour  shows  that  this  area  is 
colder  than  the  rest  of  the  system,  again  indicating  that 
heat  from  the  mixing  reactions  is  being  absorbed  by  the 
Ni  crystal  as  it  melts.  The  Ni  portion  is  completely 
disordered  at  280  ps.  The  snapshot  corresponding  to  the 
end  of  the  trajectory  (710  ps)  shows  that  the  mixing  is 
complete  and  the  Ni  and  A1  are  distributed  in  a 
homogeneous  fashion.  These  figures  indicate  that  in  this 
process,  there  are  significant  spatial  and  temperature 
inhomogeneities  in  the  material  that  will  lead  to  large, 
ever-changing  concentration  and  thermal  gradients. 

The  snapshots  in  Fig.  4  allow  an  examination  of  the 
initial  mixing  reactions  at  the  interface.  Each  snapshot  in 
Fig.  (4)  corresponds  to  three  views  of  the  same 
configuration;  the  left-most  snapshot  shows  the  bilayer 
and  the  middle  and  right-most  views  are  the  top  and 
bottom  of  the  simulation  cell.  Since  periodic  boundary 
conditions  are  imposed  in  the  simulation,  the  top  and 
bottom  show  both  sides  of  the  Ni/Al  interface, 
respectively. 
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Figure  4.  Snapshots  of  the  system  at  180  ps  from  three 
perspectives  (see  text).  The  middle  and  right-most  frames 
show  a  perspective  at  the  bilayer  interface. 

The  snapshots  correspond  to  the  material  at  180  ps, 
immediately  preceding  the  time  at  which  the  A1  melts 
(192  ps).  The  interface  contains  regions  of  both  reacted 
and  unreacted  material.  The  interface  at  192  ps  has  a 
larger  degree  of  mixing;  only  a  small  portion  of  the  A1 
layer  remains  unmixed  at  the  interface.  It  is  clear  that 
uniform  mixing  at  the  interface  does  not  occur,  and  is 
consistent  with  the  inhomogenieties  evident  in  the 
temperature  contours  of  Fig.  3. 

These  inhomogeneities  in  local  temperature  and 
concentration  gradients,  coupled  with  an  ever-increasing 
system  temperature,  will  influence  mass  and  thermal 
diffusion  effects  on  the  rate  of  reaction.  Determining 
diffusion  rates  under  these  conditions  poses  quite  a 
challenge.  However,  these  effects  are  probably  enhanced 
due  to  the  small  size  of  the  bilayer  and  separation  of  the 
two  interfaces.  We  expect  that  such  rapid  changes  in 
temperature  gradients  for  a  bilayer  having  a  significantly 
larger  separation  between  reacting  interfaces  would 
probably  not  exhibit  such  rapid  changes  in  the  thermal 
gradients.  We  are  continuing  investigations  to  determine 
the  influence  on  the  layer  thicknesses  at  various  pressures 
on  the  reaction  rates,  in  order  to  unravel  and  characterize 
the  intricate  relationship  between  pressure,  temperature 
and  system  configuration. 
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CONCLUSIONS 

We  have  begun  explorations  into  the  dynamic 
behavior  of  a  layered  Ni-Al  system  subjected  to  rapid 
thermal  loading;  we  have  demonstrated  that  pressure 
strongly  influences  the  reaction  rate,  but  that  the  overall 
mechanism  and  material  features  during  the  atomic 
mixing  and  alloying  process  do  not  differ.  This  study  has 
suggested  further  explorations  into  the  effect  on  the 
reaction  of  the  reactant  space  (i.e.  the  separation  between 
interfaces  in  a  layered  system),  as  well  as  the  dependence 
on  temperature  and  pressure  in  the  alloying  reaction. 
Information  generated  such  as  this  is  crucial  for  proper 
inclusion  in  higher  level  scales  to  adequate  capture  the 
salient  physics  and  chemistry  of  a  complex, 
interdependent  process. 
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